/* asianopt.c ---- Brittle 2009
 *
 * S0 = initial price         [5, 30]
 * v = volatility             0.3
 * r = risk free interest     0.02
 * T = mature time (in years) [0.25, 10]
 * K = strike price           [1, 100]
 * W = normal random dist.    mean=0, std_dev=1
 *
 * S[t] = S[t-1] * exp ( (r - v*v/2)*dt + v * sqrt(dt) * W )
 * where Zj is normal distributed with mean=0, std dev=1
 */

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <sys/time.h>
#include <math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

int main(void) {
  float v = 0.15f, r = 0.1f, dt = 0.01978; 
  float T = 1;
  float S0 = 100, K = 105;
  int step, mc_cnt, opt_cnt;

  float drift, vsqrdt;
  float S, S_;

  int i, j, k;
  float W;
  float path_sum;
  float call_payoff, put_payoff;
  float mc_sum, mc_call, mc_put;

  gsl_rng * rng;

  struct timeval tv1, tv2, tv3, tv4;
  long long tv_comp, tv_total;

  rng = gsl_rng_alloc(gsl_rng_mt19937);
  gsl_rng_set(rng, 0);

  scanf("%d %d", &opt_cnt, &mc_cnt);

  gettimeofday(&tv3, NULL);

  gettimeofday(&tv1, NULL);
  for (k=0; k<opt_cnt; k++) {
    scanf("%f %f %f %f %f %f", &v, &r, &S0, &K, &T, &dt);

    step = T/dt + 1;
    drift = (r - v * v / 2) * dt;
    vsqrdt = v * sqrtf(dt);

    mc_sum = 0;
    mc_call = 0;
    mc_put = 0;

    for (i=0; i<mc_cnt; i++) { // MC outer loop

      path_sum = S0;
      S_ = S0;
      for (j=0; j<step; j++) {
        W = gsl_ran_gaussian_ziggurat(rng, 1); // N(0, 1)
        S = S_ * expf ( drift + vsqrdt * W );
        path_sum = path_sum + S;
        S_ = S;
      }

      call_payoff = path_sum / step - K;
      if (call_payoff < 0) call_payoff = 0;
      mc_call = mc_call + call_payoff;
      put_payoff = K - path_sum / step;
      if (put_payoff < 0) put_payoff = 0;
      mc_put = mc_put + put_payoff;
    }

    mc_call = expf ( - r * T ) * (mc_call / mc_cnt);
    mc_put = expf ( - r * T ) * (mc_put / mc_cnt);

    printf("call pay off is %f [%08X]\n", mc_call, mc_call);
    printf("put pay off is %f [%08X]\n", mc_put, mc_put);
  }
  gettimeofday(&tv2, NULL);

  gsl_rng_free(rng);

  gettimeofday(&tv4, NULL);

  tv_comp = (tv2.tv_sec - tv1.tv_sec) * 1000000 + (tv2.tv_usec - tv1.tv_usec);
  tv_total = (tv4.tv_sec - tv3.tv_sec) * 1000000 + (tv4.tv_usec - tv3.tv_usec);
  printf("compu time = %lld usec\n", tv_comp);
  printf("total time = %lld usec\n", tv_total);

  return 0;
}
